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We present a very efficient solver for the general Anderson impurity problem. It is based on the 
perturbation around a solution obtained from exact diagonalization using a small number of bath 
sites. We formulate a perturbation theory which is valid for both weak and strong coupling and 
interpolates between these limits. Good agreement with numerically exact quantum Monte-Carlo 
results is found for a single bath site over a wide range of parameters. In particular, the Kondo 
resonance in the intermediate coupling regime is well reproduced for a single bath site and the lowest 
order correction. The method is particularly suited for low temperatures and alleviates analytical 
continuation of imaginary time data due to the absence of statistical noise compared to quantum 
Monte-Carlo impurity solvers. 

PACS numbers: 71.27.+a 



I. INTRODUCTION 

Quantum impurity models have been widely used in 
condensed matter physics. Examples comprise the study 
of the Kondo-efFect^, or of adatoms on surfaces". In par- 
ticular, the success of the dynamical mean- field theory 
(DMFT) to describe strongly correlated systems has trig- 
gered efforts to develop efficient solvers for the impurity 
problem. In DMFT, the lattice problem is mapped onto 
an effective quantum impurity problem which needs to be 
solved repeatedly to satisfy a self-consistency condition. 
Generalizations of the DMFT, such as cluster DMFT or 
dynamical cluster approximations require highly efficient 
methods for the solution of the multiorbital quantum im- 
purity problem. 

The combination of the DMFT with electronic struc- 
ture methods, such as the local density approximation 
(LDA-I-DMFT approach)'^ requires the solution of impu- 
rity models with general types of the interactions, such as 
the full Coulomb vertex. The problem has been solved us- 
ing approximate methods such as a multiorbital version 
of the fluctuation-exchange approximation (FLEX) for 
weak coupling'^. A strong-coupling solver based on an ex- 
pansion in the impurity-bath hybridization has been pro- 
posed in order to address systems with open d or f-shells 
or Mott insulators'^. Recently, next-generation Monte- 
Carlo methods provide a numerically exact solution*^'^'**. 
While being suitable for general interactions, these meth- 
ods however require a sizeable numerical effort and an- 
alytical continuation is complicated through statistical 
noise in the imaginary time data. 

In this article we propose an efficient solver for the 
impurity problem which is suitable for both weak and 
strong coupling and general interactions. It is essentially 
based on a "superperturbation" , i.e. a perturbation on 



top of a solution obtained by exact diagonalization (ED) 
for a small number of bath levels. The method alleviates 
analytical continuation due to the absence of noise and is 
suitable for the study of, e.g., multiplet effects in solids. 



II. FORMALISM 

For notational convenience, we introduce the formal- 
ism for the single-orbital case. It was introduced earlier 
in the context of lattice fermion models to include spatial 
correlations beyond DMFT in Ref. 9. A generalization 
of the underlying concepts to the multiorbital case can 
be found in Refs. 10,11. In a complementary approach, 
named dynamical vertex approximation, the single- and 
two-particle Green functions of the impurity were com- 
puted using ED^-^. The Hamiltonian of the Anderson 
impurity model (AIM) reads 



-ffAiM = ^ epbl^bpo- + ^ Vp{bl^Co- + h.c.) + i/i„t[c 
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where the local impurity degrees of freedom represented 
by c\c couple to a bath of free conduction electrons with 
dispersion via a hybridization Vp. Hmt stands for any 
local interaction. Since Haim is quadratic in the bath 
operators b\b, they can be integrated out giving rise to 
the action (henceforth we switch to the path integral rep- 
resentation): 

S[c*,c] = - ^ C*^ {{iui + H) - A^cr) Cu,a + i?int [c* , c] . 

(2) 

Here fj, is the chemical potential and the sum is over 
Matsubara frequencies ujn = {'2n+ 1)tt/(3, where /? is the 
inverse temperature. In ED the continuous dispersion 
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of the bath is approximated by a finite number of bath 
levels, which corresponds to replacing the hybridization 
function A by its discrete counterpart 



(3) 



p=i 



The problem now is to determine the parameters V^,ep 
such that A*^") is the best approximation to the origi- 
nal hybridization A. Mathematically, this corresponds 
to projecting the hybridization function onto the dis- 
crete subspacc spanned by the parameters Vp^ep. Various 
methods have been proposed for this purposc^'^'^"'. One 
way is to perform a conjugate gradient minimization of, 
e.g., the distance function 
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where the sum is over Matsubara frequencies. The pa- 
rameter k, if chosen large (e.g. fc = 3), enhances the 
importance of the lowest Matsubara frequencies in the 
minimization procedure. This parameter is particularly 
important when a small number of bath sites is used to 
approximate the original hybridization. 

The quality of the approximation can be measured by 
d and is determined by the number of bath sites used to 
represent the original hybridization. On the other hand, 
the number of bath sites is limited due to the exponential 
growth of the Hilbert space. Instead of approximating A 
by a large number of bath sites, we shall follow a different 
route and rewrite the action Eqn. 2 in terms of an exactly 
solvable part and an additional bilinear term: 

S[c*,c] = ~^c:,((zc. + M)-Ai"j)c^, + Hi„t[c*,c] 
-J2cZJAl-^ ~ A^,)c^^ . (5) 



The exactly solvable part (the first line of Eqn. 5) we will 
henceforth refer to as S^^\ Clearly, the difference = 
A^} — Ai^a- can be made arbitrarily small by including 
more and more bath sites. The main point is that the 
number of bath sites which we will employ to solve (5) 
and hence the Hilbert space is much smaller than for the 
case of conventional ED. Since S'^"-' in general is non- 
Gaussian, Wick's theorem is not directly applicable. We 
formulate a perturbation theory in D by introducing new 
fcrmionic degrees of freedom in the path integral for the 
partition function by means of the following identity: 

J exp {-ng^.Dg^,]-\f - tgZlc - c*gZlf) Hf* , f] 
= det{g^„D^^g^^y^ exp{c*D^^c) . (6) 

A similar approach has been used to derive a strong cou- 
pling expansion for the Hubbard model'^''. In the above 
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FIG. 1: Lowest order diagrams for the self-energy with 
the action given by Eqn. 9 and illustration of the Dyson 
iterations: The self-energy is obtained from by summing 
up the diagrams a)-c). From Dyson's equation, we in turn 
obtain G^ , which is subsequently used in the diagrams. 



equation g^^r = g^a = ^(cuo-c^^o-)''"'' denotes the im- 
purity Green function w.r.t. the discrete hybridization 
A("), which can be calculated exactly. With this substi- 
tution, the action becomes 

S[c\c- /*, /] = + fZ^aZlc^a + cl^gZlU 

UJ(7 

From the expression for the partition function Z = 
/ / exp(-5[c*, c; /*, /])2?[c*, c]2?[/*, /], we now in- 
tegrate out c* and c by expanding the exponential 

c,9ulf^<y) and using the fact 
that the functional integral over exp(— S**^"') produces 
correlation functions that can be obtained exactly from 
the ED, such as the two-particle Green function: 



(n) 1 / * * 

Xl234 = / CiC2C3C4exp 



(-5("))I?[c*,c] , (8) 



where we have used the shorthand notation 1 = {tJicri}. 
A compact expression for the Lehmann representation 
of the two-particle Green function is given in appendix 
A. Here we carry out this expansion up to fourth order 
(note that odd terms drop out of the expansion), with 
the result 

s^[.f\n - - E r^MVaUa + tI;! nhnu ■ (9) 



Here 7'"^ is the two-particle vertex constructed from the 
two-particle Green function as (henceforth we omit the 
superscript "(n)" on x and 7) 

71234 = 5rr333^(Xl'2'3'4' - X?'2'3'4' ).9? 2ff^3 ' (10) 
with X?234 = P{9l2gZi - 514.932). 

The auxiliary Green function is given by 

= -.9c..(.9c.. + (A£"j - A„,)-i)-i.g„, . (11) 

This function has some remarkable properties. Let us dis- 
cuss these for the case A^") =0 and Hubbard interaction 
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FIG. 2: (Color online) Imaginary part of the impurity Green 
function obtained by the superperturbation using different 
numbers of bath sites for /3 = 30 and U — 3 The representa- 
tion of the exact hybridization by A'"' is shown in the inset 
for n > (A(°> = 0). 

Hint = Un^n-i. Then we have D = —A and we can ap- 
proximate G-^ for strong hybridization by G ~ <?. In the 
weak coupling limit, i.e. U ^ 0, g approaches the bare 
Green function and 7 ~ t/, so that the dual perturbation 
theory becomes equivalent to conventional perturbation 
theory. On the other hand, for A small, itself is small 
and can obviously be approximated as G^ w guaA^a-guja- 
This generates a fast converging strong coupling pertur- 
bation expansion around the atomic limit. In fact, one 
can show that in this case the expression for Green's func- 
tion obtained by a hybridization expansion of the Green 
function given in Ref. 5 is contained in the lowest order 
diagram of our expansion. To see this, we write the aux- 
iliary Green function as G-^ ~ Gq + G^Ti^^^^G^ , where the 
self-energy correction of diagram a) in Fig. 1 is given by 
E^^^ = — 7i234(Gg)43. In order to compare this to known 
results, we have to relate the auxiliary Green function to 
the Green function of the impurity, G12 = — (ciCj). The 
fact that we have used an exact identity to introduce the 
auxiliary fermions allows us to establish an exact relation 
between G and G-^ , i.e. 

Gtjcr = + {guicrD^cr) ^ G{^„{D^aG^cr) ^ (12) 

(see Ref. 9). Inserting the approximation to G^ into 
this equation yields the following expression for G after 
some straightforward algebra: 

Gi2 = [g{Dg + l)-\2- [{I + gD)-\v 

X (x - X")i'2'3'4' [{g + D-^)-\,3, [{Dg + l)-i]2'2, 

(13) 

where we have used that fact that the vertex 7 is ex- 
pressed in terms of the two-particle Green function by 
Eqn. 10. Considering the first term in the expression 
for x° yields a contribution —Pgi2g3i[{g + -D^^)^^]43 = 
— /3Tr[g(g + Z?^^)"^]. In order to compare with Ref. 5, 
we take D — —A, as above. In this case, g corresponds to 



FIG. 3: (Color online) The contribution of different diagrams 
to the superperturbation result compared to the exact solu- 
tion (open circles), for a single bath site. The paramters are 
otherwise the same as in Fig. 2. Diagram a) yields by far the 
largest correction (upward triangles) to the initial solution 
obtained from ED (filled circles). 

the Green function for the atomic limit. For small A we 
have (1-A.g)~i 1 and (g-A^i)"^ - A. Gathering 
the results we can approximate G in the limit of small A 
as 

G12 w 312 + gi2 fi Tr[5 A] + X1234A43 , (14) 

which is essentially Eqn. (13) of Ref. 5. 

Hence the dual perturbation theory has the correct 
limiting behavior in both and weak and strong coupling 
limits. It therefore interpolates between these limits, so 
that sensible results can be anticipated even in the inter- 
mediate coupling regime. In the following we will demon- 
strate that this indeed is the case. For intermediate cou- 
pling, we further exploit the possibility to improve the 
starting point of the perturbation theory by expanding 
around the ED solution for a finite number of bath sites. 
In this case, low energy Kondo physics and the high en- 
ergy incoherent features are correctly reproduced. 

III. RESULTS 

The results shown in the following were performed us- 
ing the diagrams a) to c) depicted in Fig. 1. If not oth- 
erwise stated, the results were obtained by making use of 
the Dyson iterations illustrated in the same figure: The 
self-energy was calculated using the bare auxiliary Green 
function, Eqn. 11, on the first iteration. Inserting the 
self-energy into the Dyson equation yields a new Green 
function which is subsequently used in the diagrams on 
the next iteration. This procedure is carried out until 
self-consistency and converges in typically less than ten 
iterations. 

In the following we will discuss results obtained for the 
case of Hubbard interaction _ffi,-,t = U Jl^ dTn'j{T)n^{T). 
In order to test our approach, we performed calculations 
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for up to n = 3 bath sites. The calculations were done for 
a semielliptical density of states of bandwidth W = 4t, 
with = ~2it^/{uj + ^/uFTW). We take the half- 
bandwidth as the energy unit: W/2 = 1. 

In Fig. 2 we present results for U — 3, obtained for 
different number of bath sites up to n = 3. The quality 
of the representation of A by A*^") is shown in the inset 
for n > 0. In order to determine the parameters Vp and 
Ep, we have minimized the distance function, Eqn. 4 for 
A; = 3. This choice enhances the importance of the low- 
est frequencies in the minimization procedure. As can 
be seen in the inset, this condition results in A and A*^") 
being equal on the first n Matsubara frequencies. This 
turns out to be a good starting point for the perturbation 
theory. Using a smaller k leads to a A'-"^ which better 
represents the tail of the hybridization function and gen- 
erally leads to worse results. 

The superperturbation results are compared to those 
of a numerically exact continuous-time quantum Monte- 
Carlo (CTQMC) calculation. One can see that while the 
expansion around the atomic limit lacks accuracy, a dras- 
tic improvement occurs for a single bath site, although 
the difference between A^^^ and A is still significant. For 
n > 1, the results are essentially converged. We find qual- 
itatively the same behavior for a wide-range of U values. 
The fast convergence w.r.t the number of bath sites sig- 
nificantly reduces the computational effort compared to 
CTQMC calculations. 

In the figure we have plotted the results obtained from 
summing up skeleton diagrams. The use of skeleton di- 
agrams is theoretically relevant since in this case the re- 
sult is conserving in the Baym-Kadanoff sense^*'. The re- 
sults obtained from the first Dyson iteration or from the 
lowest-order approximation, i.e. = Gq + GqT,^Gq, 
however, achieve the same quality of approximation (not 
shown here). 

Let us now investigate the role of the different dia- 
grams in the perturbation expansion. To this end, we 
have plotted the results for the same parameters as in 
the previous figure, for different combinations of the di- 
agrams shown in Fig. 3. We have used only a single 
bath site, for which the exact solution (open circles) and 
the initial solution obtained by ED (filled circles) differ 
substantially. One can clearly see that diagram a) yields 
by far the largest correction to the initial solution. The 
value on the first Matsubara frequency is almost exactly 
reproduced. This might be expected, since A*^^^ is iden- 
tical to A on the first Matsubara frequency. The largest 
deviations occur for the second and third frequency. We 
have zoomed into this region in the inset. One can see 
that all diagrams give a correction in the right direction, 
whereby the correction by diagram c) is negligible. 

From Fig. 4 it is obvious that also spectral properties 
are correctly reproduced. We show the maximum en- 
tropy density of states^^ (DOS) obtained from the imag- 
inary time data. The analytical continuation of the quan- 
tum Monte-Carlo data (dashed-dotted line) exhibits the 
two Hubbard bands at w = ±[//2 and shows the Kondo 
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FIG. 4: (Color online) Comparison of the maximum entropy 
density of states obtained by superperturbation for no (A'-"-') 
and one bath site (A'^') to the numerically exact (continuous- 
time quantum Monte-Carlo) result. While the superpertur- 
bation around the atomic limit (A'"-*) does not reproduce the 
Kondo resonance, the perturbation around the solution for 
one bath site already contains this physics. 
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FIG. 5; (Color online) Low temperature results for ?7 = 3 and 
fi — 100. Shown are the imaginary time Green function (up- 
per and right axis) and the corresponding maximum entropy 
density of states (lower and left axis). The result obtained 
by superperturbation around two bath sites requires consid- 
erably less computational effort compared to QMC and is 
almost indistinguishable from the exact (CTQMC) result. 



resonance at the Fermi level. We cannot reproduce the 
Kondo physics by perturbing around the atomic limit 
(A'^*'\ solid line), in accordance with the findings in Ref. 
5. However, perturbation around the ED solution for 
a single bath-site already captures the Kondo resonance 
and yields good agreement compared to the exact solu- 
tion. 

In order to demonstrate that the approach also works 
for low temperatures, we present results for T = 0.01 
in Fig. 5. Although the expansion around the solution 
for a single bath site (A'^') shows small deviations in 
the imaginary time Green function G(r), the approxima- 
tion appears insufficient as seen in the density of states. 
The superperturbation around the two bath-site solution 
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however is almost exact. 

As already mentioned, the formalism of introducing 
auxiliary degrees of freedom in the path integral repre- 
sentation presented here has originally been introduced 
in the context of lattice models, termed "dual fermion 
approach"^. In that case the action is decomposed into 
an impurity part with hybridization A and a term which 
contains the bare dispersion of the lattice fermions. Since 
this decomposition is independent of the form of the hy- 
bridization, it can be represented by the discrete hy- 
bridization A^'^\ Hence the present framework can be 
used to combine dual fermion calculations with ED for 
the efficient solution of lattice models. 



IV. CONCLUSIONS 

In conclusion, we have presented an efficient approxi- 
mate solver for the quantum impurity problem. The for- 
malism straightforwardly generalizes to the multiorbital 
case and is suitable for a general form of the interaction. 
The perturbation expansion has been shown to be con- 
vergent in both weak and strong coupling limits. It was 
further demonstrated that the solver also works well in 
the intermediate coupling regime and down to low tem- 
peratures. The approximation is controlled in the sense 
that the final result can be judged on the basis of how 
well the input hybridization is approximated by a finite 
number of bath sites. The method can be used for the 
study of multiplct effects in solids in the context of real- 
istic LDA+DMFT calculations. 
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APPENDIX A: LEHMANN REPRESENTATION 
OF THE TWO-PARTICLE GREEN FUNCTION 

In this appendix, we derive a compact expression for 
the Lehmann representation of the two-particle Green 
function (2PGF). A similar expression was given in Ref.^^ 



without explicit consideration of the singular contribu- 
tions. By definition, the 2PGF in Matsubara space is 
given by 

P Jo Jo Jo 

x(r,c,(ri)4(r2)c..(r3)ct,(0)) . (Al) 

Here time translation invariance of the imaginary time 
2PGF has been used. Note that here the frequencies in 
the exponential corresponding to annihilation and cre- 
ation operators have the same sign in contrast to the 
usual definition for the Fourier transform. Correspond- 
ingly, energy conservation requires uji-\-tjj2 + u!j,-\-uji = 0. 
By restricting the range of integration such that time or- 
dering is explicit, one obtains 3! different terms. These 
can be brought into the same form by permuting the 
operators and corresponding frequencies. By the anti- 
commutation relations, each term picks up the sign of 
the permutation. After introducing the sum over eigen- 
states, the 2PGF can be written as 

ijki n 

X sgn(n)(*|OnJj) {j\OnJk) {k\OnJl) {l\cl\i) , 

(A2) 

where the first sum is over the eigenstates and the second 
over all permutations 11 of the indices {123}. We further 
have defined Oi = 0^, O2 = and (D3 = Ca' and e.g. 
Hi denotes the permutation of the first index. Here the 
different choice of convention for the Fourier transform 
simplifies the notation since otherwise the sign of the 
frequency associated with the creation operator would 
have to be permuted. The function (j) is given by the 
integral 

(f>{Ei,Ej,Ek,Ei,uji,u}2,ujs.) = 
Jo Jo Jo 

y^Q(Ek-'Ei)T3 ^ gi(wiTi +0)2 1-2+^3X3) (A3) 

The latter expression can be evaluated by taking care of 
the delta functions that arise from equal energies, with 
the final result 
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APPENDIX B: TWO- PARTICLE VERTEX IN 
THE ATOMIC LIMIT 

For the calculations without bath-sites, we have used 
an explicit expression for the two-particle vertex in the 
atomic limit. This can be obtained from result of Ap- 
pendix A together with the definition of the vertex, 
Eqn. 10. Here we reintroduce explicitly the dependence 
on W4, which yields a more symmetric form of the re- 
sult. We also switch back to the usual convention for 



the Fourier transform, for which the energy conservation 
reads oji + uj-^ ~ u)2 + oJi. 

Using that the eigenenergies oi H — fin for the impu- 
rity states |0),| Ti),! T),| i) are given by 0, 0, i7/2, i7/2, 
respectively, after some algebra the vertex is obtained as 
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